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ABSTRACT 

We investigate the final mass of the first stars by simulating the gravitational collapse of a primordial 
gas cloud until ~0.1 Myrs after the formation of the primary proto-first-star. According to our 
numerical experiment using the radiative hydrodynamics calculations, we find that the mass accretions 
onto the proto-first-stars are significantly suppressed by the radiative feedback from themselves. As a 
result, we find five stars formed in this particular simulation, and that the final mass of the stars are 
< 6OM0, including a star of AAMq. Although it is larger than 0.8 M©, presence of such a low mass 
star infer the existence of first stars in the local universe. 
Subject headings: early universe — HE regions — radiative transfer — first stars 



1. INTRODUCTION 

Formation of the first stars has been investigated in- 
tensively in the last decade mainly from theoretical as- 
pects. Following the theoretical predictions, the first 
stars form in the m ini-halos of mass ^ 10^ — lO^Mc?) 
(iHaiman et al.l 119961: F Tegmark et all 119971: iNishi fc Susal 
19991: iFuller fc Co uchman 2000' 'Abel Brvan, Norman' 
Bromm. C oppi & Larson 2002; Yoshida et al. 
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The ingredient of first star formation is the primordial 
gas, which does not contain heavy elements or cosmic 
dusts. Because of the lack of these efficient coolants, the 
primordial gas cools inefficiently especially at low tem- 
peratures (T < lO^K). Therefore, the gas is kept rela- 
tively warm (~ lO^K) while it collapses to form stars, in 
contrast to the case of the interstellar gas, whose tem- 
perature is ~ lOK during the present-day star forma- 
tion for nn ^ lO^^cm"^ (e.g. Omukai 2000). As a re- 
sult, the gravitationally collapsing primordial clouds are 
very massive ~ lOOOM©, since they have to be more 
massive than the Jean's mass, which is proportional to 
T^/^. In addition, formation of such a massive prestellar 
core leads to huge mass accretion rate onto the protostar 
in the later mass accretion phase f e.g. Omukai & Nishi 
119981 : lAbel. Brvan. NormanI [200I l^hida et al. 2006). 
Following these theoretical evidences, the first stars were 
once expected to be very massive (> lOOM©). 

On the other hand, the studies on the mass accretion 
phase have advanced recently, which revealed that the 
first star formation pr ocess seems to be m ore complicated 
than expected be fore (Cl ark et ahl l2011al lb: Smith et al.l 
[20TTI : IGreif et al.l[201L .20121 ). They found that a heavy 
disk formed around the primary protostar, because of 
the angular momentum of the prestellar core, gained by 
the tidal interactions with other cosmological overdense 
regions. The heavy disk fragments into small pieces since 
it is gravitationally unstable. As a result, a "star cluster" 
could be formed instead of a single very massive star. 

These results seem to be robust until the primary pro- 
tostar grows to > 2OM0. After the mass of the pro- 
tostar exceeds ~ 2OM0, significant ultraviolet radiation 
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flux will be emitted from the protostar (jOmukai fc Pal 



120031 : iHosokawa fc Omukaill2QQ9l : iHosokawa et al.ll2011[ ) 
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Thus, the ultraviolet radiation from the protostar signif- 
icantly affects the later evolution of the system. 

Hosokawa et al. (2011) have addressed this feedback 
effect directly by two dimensional radiation hydrody- 
namics simulations. They found that the accretion disk 
around the primary protostar is photoevaporated due to 
the radiative feedback, followed by the rapid decline of 
the mass accretion rate onto the protostar. The flnal 
mass of the protostar in their simulation is 43 Mq. 

Stacy et al. (2012) also tried to assess the flnal mass of 
the flrst stars in their three dimensional cosmological cal- 
culations including the fragmentation of the disk as well. 
They also found that the ultraviolet radiative feedback 
strongly suppress the mass accretion onto the protostars. 
However, the integrated physical time is ~ SOOOyrs, 
which is too short to predict the flnal mass of the flrst 
stars, since it will take ~ O.lMyr s until the protostar 
settles onto the main sequence (e.g. IHosokawa fc Omukail 
l2009|). 

In this paper, we report the results of three dimensional 
radiation hydrodynamics simulations on the formation of 
flrst stars that follow the evolution of the system for 0.1 
Myrs after the formation of the primary protostar. We 
take into consideration the three dimensional effects, as 
well as the radiative feedback from the protostars. 

2. NUMERICAL SIMULATION 

We study the formation of flrst stars by radiation hy- 
drodynamics simulations. We employ the Bonner-Ebert 
sphere of nu = 10^cm~^,T = 200K as the initial condi- 
tion of the simulation. This initial condition is mo tivated 
by t h e cosmological simul ations (j Abel. Bryan. NormanI 
I2OO2I : I Yoshida et al. I l2003f ) in which such clouds are 
found in minihalos of mass 10^ — lO^M©, at the "loi- 
tering" phase. The loitering phase corresponds to the 
epoch when the cloud becomes quasi-static, because H2 
line cooling for nn ^ lO'^cm"^ becomes less efflcient than 
that for Tin < lO^cm"^. In order to make the gas cloud 
slightly gravitationally unstable, we increase the den- 
sity by 20%. As a result, the total mass of the cloud 
is 26OOM0. We also add uniform rotation to the gas 
sphere with ^^ = 2 x 10~^^rad/s. The angular velocity 



of this value results in very similar specific angular mo- 
mentum distri bution to that found in the cosmological 
simulations by lYoshida et al.l ()2QQ6f ) at the final stage of 
the run- way collapse phase. 

We use the c ode Radiation-SPH (Susa fc Umemural 
I2QQ4I : iSusal l2QQ6l ) in order to solve the equations of hy- 
drodynamics, non-equilibrium primordial chemistry of 
six species, e~, H+, H, H~, H2", and radiation trans- 
fer of ultraviolet photons. We solve the propagation of 
Lyman- Werner (LW) band photons as well as the ioniz- 
ing photons. We use updated self-shielding function for 
LW radiation transfer (Wolcott- Green et al.i l2Qllf ). The 
mass of an SPH particle in the present work is set to be 
^^SPH = 4.96 X 10~^Mq and the number of neighbor par- 
ticles is A^neib = 50, that correspond to the mass resolu- 
tion of Mres = 2A^neibrnspH = Q.496M0 (lBate fc BurkertI 
[1997). This mass resolution is equivalent to the Jeans 
mass o f nn = lQ^^cm~^, T = 300K, and comparable to 
that in lStacvet all (j2Q12f ). 

In order to trace the evolution in mass accretion phase, 
we employ sink particles. If the density at an SPH par- 
ticle exceeds ngink = 3 x lO^^cm"^, we change the SPH 
particle into a sink particle. In addition, if SPH parti- 
cles fall within the sphere with radius of race = 30AU 
centered on a sink particle, and they are gravitationally 
bound with each other, these SPH particles are merged 
to the sink particle, conserving the linear momentum and 
mass. The accretion radius rac e is a gain comparable to 
that employed in IStacv et al.l (j2Q12l ). We remark that 
the sink-sink merging is not allowed in the present nu- 
merical experime nt, since the radius of the p rotostar is 
less than ~ lAl KjHosokawa fc Omukalll2QQ9f ). which is 
much smaller than the employed accretion radius race- 
We regard the mass of the sink particles as the mass of 
protostars. We also assume that the sink particles do 
not push surrounding SPH particles, i.e. the pressure 
forces from sink particles to surrounding SPH particles 
are omitted. The recipe of sink particles employed in the 
present work is known to overestimat e the ma ss accre- 
tion rate(e^£^_^a^e^_al. 1995; Bromm. Coppi fc Larson 
I2QQ2I : iMartel et aDl2QQ6D . In addition, the employed ac- 
cretion radius race is 30 AU, which is much larger than the 
radius of protostars (jHosokawa fc Omukail |2QQ9[ ) . Thus, 
we have to keep in mind that resultant mass of the formed 
sink particles would be larger than the actual mass of first 
stars. We also remark that we cut the central spherical 
region with radius 0.6pc out of the cloud, just after the 
formation of the first sink in order to save the compu- 
tational time. The outer envelope of r > 0.6pc hardly 
affect the inner region within lO^yrs. 

The luminosity and the effective temperature of the 
protostar is obtained by int erpolating the results from 
iHosokawa fc Omukail (j2QQ9l ). which are tabulated as 
functions of protostellar mass (M*) and the mass accre- 
tion rate (M*). 

3. RESULTS 

We perform radiative hydrodynamics simulations of 
the first star formation with radiative feedback. We also 
perform a run with no feedback for comparison. 

3.1. Fragmentation of the disk around the primary 
protostar 



We start the simulation from the rigidly rotating 
Bonner-Ebert sphere around the loitering phase. The 
cloud starts to collapse in run-away fashion, i.e. the cen- 
tral density keeps growing while the outer part of the 
cloud is left in the envelope. As a result, density pro- 
file of (X r~^-^ is b uilt up during the run- away phase 
(jYoshidaet al.ll2QQ6f ). 

Eventually, the central density exceeds ngink, and a 
sink particle is formed at the center of the cloud. The 
surrounding gas starts to accrete onto the sink parti- 
cle subsequently. Since the gas has a significant amount 
of specific angular momentum, the accreting gas forms 
an accretion disk around the sink particle. The amount 
of specific angular momentum in the run-away phase is 
close to that of the similarity solution, which is approxi- 
mately 0.5 times the value of Kepler rotation at th e Jeans 
radius, i.e. the core radius (Yoshida et al."2006). Thus, 
the radius of the disk is 0.25 times smaller than their orig- 
inal radius at the run- away phase, since the centrifugal 
force is proportional to the square of the specific angular 
momentum. 

After the formation of the small gas disk around the 
first sink, the gas keeps accreting onto the disk. As a 
result, the mass and the radius of the disk increase. At 
the same time, the temperature of the disk decreases by 
the radiative cooling. The left column of Fig(T] shows 
the edge-on views of the disk column density at three 
epochs corresponding to 320yrs, 620yrs and 860yrs after 
the formation of the first sink. The red crosses denote 
the position of sink particles. 

In the early epoch of the accretion phase, a smooth 
disk forms around the sink particle (top panel), followed 
by the formation of the spiral arms (middle panel) and 
the fragmentation of the arms (bottom panel). We can 
also find a few high column density peaks in the bottom 
panel. In fact, next sink particles are born from these 
peaks within a few hundred years. 

The right column of Fig(T] shows the color contour of 
Toomre's Q-parameter, which is given as 

Q = Cs^orh/{7TGa) 

in case we assume Keplerian motion. Here Cg denotes the 
sound velocity of gas, l^orb is the angular velocity around 
the central sink particle and a is the column density of 
the disk. In the case that the Q-parameter is less than 
unity, a smooth disk with density perturbations becomes 
gravitationally unstable. In fact, the Q-parameter of the 
disk at the early phase (top) is already less than unity, 
so the disk is unstable (middle and bottom). 

The time scale of disk f ragmentation is given by the 
linear perturbation theory ()Toomrelll964f ) , which is given 
as ^OTb(Q~^ — 1)~^/^. The typical value of Q-parameter 
in the disk at early phase (top) is ~ 0.5, and the angu- 
lar velocity is florh — 10~^s~^. Thus the perturbation 
growth time scale is ~ 20yrs, which is comparable to the 
time scale of the generation of the spiral structure. Thus, 
spiral structures develop due to the gravitational insta- 
bility, which could be understood as Toomre's criteria. 
On the other hand, it takes several hundred years for an- 
other sink to be born (bottom), and seemingly via the 
fragmentation of the spiral arms. Therefore, formation 
of the sink particles in the disk cannot be understood 
solely by the simple Q-parameter argument above, but 
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Fig. 1. — Edge-on view at the begging of the accretion phase. Left column: Three snapshots (320yr, 620yr, 860yr from top to bottom) 
of color contour of the gas column density. Right column: Toomre's Q-parameter at the same moment. 



requires non-linear calculations. 

3.2. Effects of radiative feedback 

Fig 12] illustrates the evolution of the edge-on view of the 
central 10^ AU (O.OSpc) in radius. The color shows the 
fraction of H2 molecules, and the transparency denotes 
the gas density. Small spheres are the position of sink 
particles. Initially, H2 fraction is quite high {yn^ ~ 10~^, 
upper left panel). Eventually, the polar region is pho- 
todissociated as the sink particles grow (top right). Some 
H2 rich regions remain along the equatorial plane due to 
the self-shielding (bottom left), but finally they disap- 
pear after O.lMyrs (bottom right). 

Fig |3] illustrates the evolution of the system on density- 
temperature plane. Color map shows the frequency dis- 
tribution of SPH particles on the plane. Four panels show 
the snapshots at t =-10yr, 2450yr, 5510yr, 100250yr, 
respectively. The distribution of SPH particles on the 
top left panel is similar to the well known curve of the 
evol ution of collapsing primordial gas in run- away phase 
(e.g. iPalla et al]ll983l ). since it corresponds to the epoch 
just before the first sink formation. After the first sink 
formation, dense gas (> lO^^cm"^) is splitted to high 
temprerature gas (< 7000K) and low temperature gas 
(< lOOOK, top right panel). The former corresponds 
to the radiatively heated gas and the shock heated gas, 
while the latter is the self-shielded cold gas orbiting 
around the sink particles. Then, the high temperature 
gas in the high density region expands due to the in- 
creased thermal pressure, generating a shock propagat- 
ing to low density region (bottom left panel). In fact, we 
also have seen the shock wave in the bottom left panel of 
FiglSl marked by white dashed curve. Finally, the dense 
cold gas disappears (bottom right), which means that 
star formation no longer proceeds. 

The solid line in Fig J3.2l shows the time evolution of 
the total mass in the sink particles in the feedback run, 
whereas the dashed line represents that without radiative 
feedback effects. It is clear that the mass accretion onto 
the sink particles is highly suppressed by the radiative 
feedback. The total sink mass of the run with feedback 
at the end of the simulation ('^0.1 Myr) is less than a 
third of that without feedback. Fig H] illustrates the time 
evolution of the mass of the each sink particle. The red 
curves represent those in the run with the feedback, while 
the green curves are those results without feedback. In 
the feedback run, we have one star more massive than 
50 Mq at 0.1 Myr (~ STM©), whereas we have three stars 
in the range of 50 — 2OOM0 in the no feedback run. 

The minimal sink mass in the feedback run is 4.4M0, 
and the total number of the sink particles is 5. On the 
other hand, in the no feedback run, the minimal sink 
mass is O.84M0, and the total number of sink particles 
is 10. It is also worth noting that fragmentation of the 
disk in the no feedback case continues until much later 
time (~ 4 X lO^yr) than that in the feedback run (< 
ISOOyr), because the molecular rich disk is not destoried 
in the no feedback run. The difference in the number 
of sink particles and the minimal mass might come from 
such effect. However, it is premature to draw definitive 
conclusion on this issue, since our results are based upon 
only a single realization. 

3.3. Es capers 



In the simulation with feedback, we find a sink particle 
is kicked away from the central dense region via the gravi- 
tational N-body interaction, so-called "slingshot" mecha- 
nism. Figinishows the trajectories of all the sink particles 
within (2 x 10^ AU)^ box at the central region. It is clear 
that one escaping sink goes away from the central region 
(to the bottom of the panel), whereas the others remain 
within the box. Such a phenomenon has already been 
reported b y previous groups without radiative feedback 
effects(e.g. Smith e t al.|[2QlT ). Thus, we confirm the the- 
oretical existence of such escapers also in our numerical 
model with radiative feedback. The velocity of this es- 
caping sink is ~ 4km/s at 0.1 pc distant from the center 
of the cloud, which is merginal to escape from the host 
minihalo of lO^M©. 

The mass accretion onto this escaping particle almost 
stops after the ejection. Consequently, the mass of the 
sink is 4.4M0, that are much smaller than the "conven- 
tional" first stars of mass > lOOM©. It is also woth 
noting that the orbit of the other sinks are excited by 
the N-body interaction with each other, although they 
are not ejected. Thus, some of the sinks going through 
relatively low density regions, which results in low mass 
accretion rate onto these sinks. 

We also remark that a star less massive than O.SMq 
is found in our higher resolution run (M^es = O.IMq) 
with feedback, although the integrated physical time 
is ~ 2 X lO^yrs (Umemura et al. 2012). Considering 
the fact that the mass resolution and the accretion 
radius of the present simulations are ~ O.SM© and 
30AU, and other higher resolution studies with/without 
feedback effects report the ejection of even lower 
ma ss stars (Umemura et al. 2012; Clark et al. 2011a. b; 
ISmithet al.ll201ll : iGreif et al.ll201lL l2012h . we presume 
that low mass stars less massive than O.SM© are born 
among the first stars, and survive through the entire his- 
tory of the universe. 

4. DISCUSSIONS 

In the presence of radiative feedback, the gas in the 
neighborhood of protostars is heated up significantly. 
This heating process occurs mainly through photodis- 
sociation of H2 molecules: Formation of H2 molecules 
works as a heating process of the gas, since formation 
process such as 3H^H24-H or H~+H ^ H2 + e~ re- 
leases the latent heat. In the absence of radiative disso- 
ciation process, collisional dissociation processes absorb 
thermal energy, which balance with the formation heat- 
ing. Thus, the net increase of H2 molecules causes ef- 
fective heating of the gas, while the net decrease of H2 
results in cooling. On the other hand, in the presence of 
strong photodissociative radiation, it overwhelms other 
collisional dissociation processes. The photodissociation 
process do not absorb thermal energy, since the energy 
required to dissociate H2 molecules is supplied by the 
radiation. Therefore, H2 dissociation is no longer a cool- 
ing process. As a result, H2 formation heating proceeds 
without hinderance in the absence of the counter process. 

In fact. Fig l6] illustrates the ratio between the H2 for- 
mation heating rate and the adiabatic heating rate as 
functions of gas temperature and density at t = 2450yr. 
It is clear that H2 formation heating is the dominant 




Fig. 2. — Edge-on views of gas distribution inside r < IC^AU (O.OSpc) at four snapshots. Left column: from top to bottom, t =- 
10yr,1160yr. Right column:^ = 5120yr, 100250yr. t represents the time after the first sink formation. Color shows the H2 fraction, and the 
small spheres represent the positions of sink particles. White arrow and dashed curve in the top right panel denotes the position of the 
shock front. 
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Fig. 3. — Four snapshots on the density-temperature plane. Color 
represents the number of SPH particles drop on the logarithmic 
bin on the plane. Thick solid lines show the resolution limit of this 
simulation, while dashed lines represent the number density above 
which the SPH particles are converted to sink particles. 
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Fig. 4. — The mass of individual sink particles are plotted as 
functions of time after the first sink formation. Red curves for the 
runs with feedback, while the green curves for the case without 
feedback. 
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Fig. 5. — The trajectories of sink particles in edge-on view. 

heating process at high densties (nn ^ lO^^cm"^, right 
panel). The temperature of these chemo- heated high 
density regions is mainly around lOOOK, and the chem- 
ical heating is also important for high temperature re- 
gion at 3000K < T < 7000K (left panel). It is also 
clear that chemical heating rate is not negligible com- 
pared to the adiabatic heating rate even at lower densities 
(~ lO^cm"^) and lower temperatures (T ~ 300K), which 
corresponds to the dissociated polar regions. Thus, 
chemical heating of H2 formation play important roles 
on the dynamical evolution of this system. 

Heating through photoionizati on is also important. I n 
fact, similar calculation in 2D bv lHosokawa et al] ()2Qlll ). 
photoionization is the dominant heating process. They 
found that gas in the polar region is highly ionized and 



heated up to > lO^K. Finally, the disk is photoevapo- 
rated mainly due to the photoheating process through 
ionization. In the present calculation, however, ioniza- 
tion is not the dominant heating process. The reason 
simply comes from the fact that the spatial resolution 
of the present simulation is not e nough to resolve the 
"initial" Stromgren sphere (jSpitzerl ll978) . The SPH par- 
ticles in the very vicinity of the source star is heated by 
ionization only "mildly", T < lO^K, which cannot cause 
the drastic propagation of D-type ionization front which 
is found in Hosokawa et al. (2011). In other words, this 
radiative feedback effect is underestimated in the present 
simulation. 

As noticed in section [21 the algorithm of mass accre- 
tion onto sink particles overestimates the mass accretion 
rate. Combined with the fact that feedback is underes- 
timated, the mass accretion rate obtained in the present 
numerical experiment is overestimated at any hand. Ac- 
cordingly, the final mass of the sinks should be regarded 
as an upper limit of the actual mass of first stars. If we 
could perform simulations at higher resolution with more 
realistic mass accretion conditions, the final mass could 
be smaller than ~ SOM©. We also should keep in mind 
that this result comes from a numerical experiment of one 
realization. Thus, this upper limit should be regarded as 
a guide, since a slightly different initial condition could 
cause different result because of the chaotic nature of the 
system. In any case, the final mass of the sinks seem not 
to exceed IOOMq from typical initial conditions. On the 
other hand, at least one of the protostars cannot be less 
massive than M ~ 2OM0, since the radiative feedback 
becomes prominent only for M* > 20Mq. Thus, the 
mass of the primary star among the first stars formed in 
a single mini-halo will fall in the range of 20 — 60Mq. 

The initial condition of the present experiment is a 
rigidly rotating Bonner-Ebert sphere. Although its an- 
gular momentum distribution just before the sink for- 
mation is close to that of cosmological simulations, it is 
not fully cosmological. In cosmological simulations, the 
direction of angular momentum of the disk around the 
primary protostar changes depending on the stages, since 
the gas motion is more turbulent. In addition, we need 
more numerical experiments starting from various initial 
conditions, in order to obtain statistical quantities such 
as the initial mass function (IMF) of first stars. Thus it 
is important to perform fully cosmological simulations of 
this sort, which is left for future works. 

5. SUMMARY 

In this paper, we investigated the suppression of mass 
accretion onto proto-first-stars by the radiative feedback 
from themselves. We performed numerical experiment 
of the formation of first stars using three dimensional ra- 
diative hydrodynamics code RSPH combined with sink 
particle technique. Consequently, we find that the mass 
accretion is suppressed significantly and the final out- 
come is a multiple stellar system consisting of five stars 
of 1 — 60Mq. The fact that low mass stars are found in 
this work infer the possible existence of first stars in the 
local universe, although the mass of the formed stars in 
our simulation is still larger than 0.8 Mq. 

We thank T. Hosokawa for providing the data of proto- 
first-stars. We also thank K. Omukai for careful reading 
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Fig. 6. — The ratio between the chemical heating rate and the adiabatic heating rate is plotted as functions of temperature (left panel) 
and density (right panel) at t = 2450yr. 
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